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Abstract. This article presents a set of spectro-goniometric 
measurements of different water ice samples and the compar¬ 
ison with an approximated radiative transfer model. The ex¬ 
periments were done using the spectro-radiogoniometer de¬ 
scribed in Brissaud et al.| ( |2004| ). The radiative transfer model 
assumes an isotropization of the flux after the second inter¬ 
face and is fully described in |Andrieu et al.|p015| |. 

Two kind of experiments were conducted. First, the spec¬ 
ular spot was closely investigated, at high angular resolution, 
at the wavelength of 1.5 pm , where ice behaves as a very 
absorbing media. Second, the bidirectional reflectance was 
sampled at various geometries, including low phase angles 
on 61 wavelengths ranging from 0.8 pm to 2.0 pm . 

In order to validate the model, we made a qualitative test 
to demonstrate the relative isotropization of the flux. We also 
conducted quantitative assessments by using a bayesian in¬ 
version method in order to estimate the parameters (e.g. sam¬ 
ple thickness, surface roughness) from the radiative measure¬ 
ments only. A simple comparison between the retrieved pa¬ 
rameters and the direct independent measurements allowed 
us to validate the model. 

We developed an innovative bayesian inversion approach 
to quantitatively estimate the uncertainties on the parameters 
avoiding the usual slow Monte Carlo approach. First we built 
lookup tables, and then searched the best fits and calculated 
a posteriori density probability functions. The results show 
that the model is able to reproduce the spectral behavior of 
water ice slabs, as well as the specular spot. In addition, the 
different parameters of the model are compatible with inde¬ 
pendent measurements. 


1 Introduction 


Various species of ices are present throughout the Solar Sys¬ 
tem. From water ice and snow on Earth to nitrogen ice on 
Triton ( |Zent et ah} |1989| , not forgetting carbon dioxide ice 
on Mars ( (Leighton and Murray) |1966| . Ice an snow cov¬ 
ered areas have a strong impact on planetary climate dy¬ 
namics, as they can lead to significant regional scale albedo 
changes at the surface and surface/atmosphere volatiles in¬ 
teractions. The physical properties of the cover have also an 
impact on the energy balance : for example the albedo de¬ 
pend on the grain size of the snow ( [Dozier et al.| [3009 ; Negi 
|and K okhanovsky, 2011), on the roughness of the interface 
( Fhermitte et al., 2014) on the presence or not and the phys¬ 
ical properties of impurities ( Dumont et al.| 20l4|), or on the 
specific surface area ( [Picard et al.| 2009; Mar y et al.||2013| ). 
The study and monitoring of theses parameters is a key to 
constrain the energy balance of a planet. 

Radiative transfer models have proven essential to retrieve 
such properties and their evolution at a large scale, and dif¬ 
ferent families exist. Ray tracing algorithms, such as those 
described in |Picard et~ah| ( |2009| ) for snow or Pilorget et al. 


( 2013| ) for compact poly crystalline ice simulate the com 
plex path of millions of rays into the surface. They pro¬ 
vide very accurate simulations, but have the weakness of 
being time consuming. Analytical solutions of the radiative 
transfer in homogeneous granular media have been devel¬ 
oped for example by Shkurato v et al.| ( |1999| ) and |Hapke 
( |1981| ). They are fast, and when the surface cannot be de¬ 
scribed as homogeneous, they must be combined with an¬ 
other family of techniques such as discrete ordinate meth¬ 
ods like DISORT ( [Stamnes et al.[|1988] ). These methods have 
been widely studied on Earth sno w ([Carmagnola et al.| 2013 
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2 


Andrieu: Radiative transfer model for contaminated slabs : experimental validations 


-2d 

*+ -Ice matrix 


h 


Granular 
subtrate 

Figure 1. Scheme of the surface representation in the radiative 
transfer model applied on the laboratory measurements, h repre¬ 
sents the slab thickness. 0 represent the mean slope to describe the 
surface roughness. 
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2004) and other planetary cryospheres ( Appe re et al.||20ll] 


Eluszkie wicz and Moncet| |2003| ), modeling a granular sur¬ 


face. Compact poly crystalline ices have however been rec¬ 


ognized to exist on several objects : CO 2 on Mars (Kieffer 


land Titus[|2001[ |Eluszkiewi cz et al.[|2005| ), N 2 on Triton and 
Pluto, ( |Zent et ah 1989| |Eluszkiew icz and Moncet| |2003| ) 
and probably SO 2 on Io (Eluszkiewicz and Moncet 2003), 
owing to the very long light path-lengths measured, over sev¬ 
eral decimeters. Radiative transfer in compact slabs is differ¬ 
ent from in granular media. We developed an approximated 
model ( [Andrieu et al.||2015] ) designed to study contaminated 
ice slabs, with a fast numerical implementation, that has al¬ 
ready been numerically validated. The main objective of the 
model is the analysis of massive spectro-imaging planetary 
data of these surfaces. For this purpose, it is semi analytic 
and fast implemented. It is designed to retrieve the variations 
of thickness and impurity content of compact polycrystalline 
planetary ices. 

In the present article, we will test the accuracy of this 
approximated model on laboratory spectroscopic measure¬ 
ments of pure water ice Bidirectional Reflectance Distribu¬ 
tion Function (BRDF). The goal is to propose an inversion 
framework to retrieve surface properties, including uncer¬ 
tainties in order to demonstrate the validity of the approach. 
In order to speed up the inversion, we based the algorithm on 
look-up tables that minimize the computation time of the di¬ 
rect model. This strategy will be very useful to analyze real 
hyperspectral images. The thickness of ice estimated from 
the inversion is validated in comparison to real direct mea¬ 
surements. Also the specular lobe is adjusted to demonstrate 
that the model is able to reasonably fit the data with a coher¬ 
ent roughness value. 


2 Description of the model 


The model, Andrieu et al 


IpoBI , isi 
>ke (19811) a: 


inspired from an exist¬ 


ing one described in Hapke (1981) and Doute and Schmitt 



substrate 


Figure 2. Illustration of the radiative transfer in the surface. 
Anisotropic transits are represented in red. F is the incident radia¬ 
tion flux, Rspec and RDiff are respectively the specular and diffuse 
contributions to the reflectance of the surface. r s is the lambertian 
reflectance of the granular substrate, and Ro and To are total reflex¬ 
ion and transmission factors of the slab layer. A prime indicates an 
anisotropic transit. The reflection and transmission factors are dif¬ 
ferent in isotropic or anisotropic conditions. The granular and slab 
layers are artificially separated in this figure to help the understand¬ 
ing of the coupling between the two layers. On the top : illustration 
of the reflections and transmission at the first interface, used in the 
calculations of R sp ec and the determination of the amount of energy 
injected in the surface, z is the normal to the surface, W / the local 
normal to a facet, i and e are respectively the incidence and emer¬ 
gence angle, and e/ is the local emergence angle for a facet. Each 
different orientation of a facet will lead to a different transit length 
in the slab. A more detailed description can be found in |Andrieu| 
|etaLl ( [20T5l ). 


( |1998| ), that simulates the bidirectional reflectance of strat¬ 
ified granular media. It has been adapted to compact slab, 
contaminated with pseudo-spherical inclusions, and a rough 
top interface. In the context of this work, we suppose a layer 
of pure slab ice, overlaying an optically thick layer of granu¬ 
lar ice, as described on figure[T] The roughness of the first in¬ 
terface is described using the probability density function of 
orientations of slopes defined in |Hapke| ( |1984| ). This distribu¬ 
tion of orientations is fully described by a mean slope param¬ 
eter 9. The ice matrix is described using its optical constants 
and it thickness. Within the slab, the model can also incorpo¬ 
rate inclusions, supposed to be close to spherical and homo¬ 
geneously distributed inside the matrix. They are described 
by their optical constants, their volumetric proportions and 
their characteristic grain-sizes. There can be several different 
types of inclusions. Each type can be of any material, except 
the one constituting the matrix : it can be any other kind of 
ice, mineral, or even bubbles. 
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Figure [2]illustrates the general principle of the model. The 
simulated bidirectional reflectance results from two separate 
contributions: specular and diffuse. The specular contribu¬ 
tion of a measurement is estimated from the roughness pa¬ 
rameter, the optical constants of the matrix, and the apertures 
of the light source and the detector. The surface is consid¬ 
ered as constituted of many unresolved facets, which orien¬ 
tations follow the defined probability density function. The 
specular reflectance is obtained integrating every reflection 
on the different facets. The total reflection coefficient at the 
first rough interface is obtained by integrating specular con¬ 
tributions in every emergent direction, at a given incidence. 
This gives the total amount of energy transmitted into the sys¬ 
tem constituted of the contaminated slab and the substrate. 
The diffuse contribution is then estimated solving the radia¬ 
tive transfer equation inside this system under various hy¬ 
pothesis. It is considered that: (i) the first transit through the 
slab is anisotropic due to the collimated radiation from the 
source, and that there is an isotropization at the second rough 
interface (i.e. when the radiation reach the semi-infinite sub¬ 
strate). For the refraction and the internal reflection, every 
following transit is considered isotropic, (ii) the geometrical 
optics are valid. This means that the size of the inclusions 
and the thickness of the slab layer must be larger than the 
considered wavelength, (iii) the inclusions inside the matrix 
are close to spherical and homogeneously distributed. The re¬ 
flexion and transmission factors of the layers are obtained us¬ 
ing an analytical estimation of Fresnel coefficients described 
in |Chandrasekhar|fl960| ); [Doute and Schmitt ( |1998 ), and a 


Andrieu et al. 


pOT5| . 


simple statistical approach, detailed in 
The contribution of the semi infinite substrate is described by 
its single scattering albedo. Finally, as the slab layer is under 
a collimated radiation from the light source, and under a dif¬ 
fuse radiation from the granular substrate, the resulting total 
bidirectional reflectance is computed using adding doubling 
formulas. 


3 Data 


3.1 Spectro-radiogoniometer 


The bidirectional reflectance spectra were measured using 
the spectrogonio-radiometer from IPAG fully described in 


Brissaud et al. ([2004]). We collected spectra in the near in¬ 
frared at incidences ranging from 40° to 60°, emergence an¬ 
gles from 0° to 50°, and azimuth angles from 0° to 180°. 
The sample is illuminated with a large monochromatic beam 
(divergence < 1°) an the near infrared spectrum covering 
the range from 0.800 pm to 4800 pm is measured by an 
InSb photovoltaic detector. This detector has a nominal aper¬ 
ture of 4.2°, that result in a field of view on the sample of 
approximately 20 mm in diameter. The minimum angular 
sampling of illumination and observation directions is 0.1°, 
with a reproducibility of 0.002°. In order to avoid azimuthal 
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Figure 3. (a) Reflectance factor at a wavelength A = 1.4pmversus 
phase angle for snow only (black crosses) and the same snow, but 
covered with a 1.42±0.27mm water ice slab (red squares). The thin 
layer of slab ice not only lower the level of reflectance as expected, 
but also seems to isotropize the reflected radiation. It is clearer on 
plot (b), that represents the same data, but normalized by the value 
at a phase angle a = 20°. This data shows that even a very thin 
layer of ice has a strong effect on the directivity of the surface. This 
justify the approximation of isotropization at the second interface 
supposed by the model, and the description of the substrate using 
only its single scattering albedo. 


anisotropy, the sample is rotated during the acquisition. The 
sample rotation axis may be very slightly misadjusted, result¬ 
ing on a notable angular drift on the emergence measured up 
to 1°. 
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3.2 Ice BRDF measurements 

The ice samples were obtained by sawing artificial columnar 
water ice in sections. These sections were put on top of an 
optically thick layer of compacted snow, that was collected 
in Arselle, in the French Alps. The spectral measurements 
were conducted in a cold chamber at 263K. Still, the ice and 
the snow were unstable in the measurement’s environment, 
due to the dryness of the chamber’s atmosphere. The grain- 
size of the snow showed an evolution and the thickness of a 
given slab showed a decrease of 0.343mm/day. Each sample 
needs a acquisition time of 10 hours. For each measurement, 
the ice slab was sliced, and its thickness was measured in 
five different locations. It was then set on top of the snow 
sample, and this system was put into the spectro-goniometer, 
and put into rotation for the measurement. As the surface is 
not perfectly plane, the measured thickness is not constant. 
This results in an 2 a standard-deviation in the measurement 
of the thickness than ranges from 0.54mm to 2.7mm in our 
study, depending on the sample. 

3.2.1 Specular contribution 

The specular contribution was measured on a 12.51 mm thick 
slab sample on top of Arselle snow. This sample is the same 
sample that sample 3 described in the next paragraph. The il¬ 
lumination was at an incidence angle of 50°, and 63 different 
emergent geometries were sampled, ranging from 45° to 55° 
in emergence, and from 170° to 180° in azimuth. A measure 
at the wavelength of 1.5 /im is represented on figure |5] (a). 
The sampling is 1° in emergence and azimuth within 47° and 
53° in emergence and 175° and 180° in azimuth. 

3.2.2 Diffuse reflectance spectra 

The diffuse contribution was measured on three samples 
of different slab thickness. The three thicknesses were 
measured on different locations of the samples with a 
caliper before the spectrogonio measurement, resulting in 
hi = 1.42 ±0.47mm, h 2 = 7.45 ± 0.84mm, h 3 = 12.51 ± 
2.7mm, respectively for samples 1, 2 and 3, with er¬ 
rors at 2a. 61 wavelengths were sampled ranging from 
0.8 /im to 2.0 /im. Spectra were collected on 39 different 
points of the BRDF, respectively for the incidence, emer¬ 
gence and azimuth angles : [40°,50°,60°], [0°,10°,20°] and 
[0°,45°,90°, 140°, 160°, 180°]. This set of angles results in 
only 39 different geometries because the azimutal angle is 
not defined for a nadir emergence. 

3.2.3 Snow diffuse reflectance spectra 

Diffuse reflectance spectra of natural snow only were also 
measured. The objective was to estimate the effect of a slab 
layer on the BRDF. Figure [3] show on the same plots the re¬ 
flectance factor versus phase angle of the snow and the snow 
covered with a 1.42 mm thick ice slab (sample 1). It illus¬ 


trates the two most notable effects of a thin layer of slab ice 
on top of an optically thick layer of snow. The most intuitive 
effect is to lower the level of reflectance : it is due to ab¬ 
sorption during the long optical path lengths in the compact 
ice matrix. The second is that the radiation is more lamber- 
tian than the one of only snow. This data gives credit to the 
first hypothesis of isotropization of the radiation formulated 
in the model (see Section [2]). The description of the bottom 
granular layer as isotropic, defined only its single scattering 
albedo, may be considered brutal, but this dataset shows that 
a thin coverage of slab ice even on a very directive material 
such as snow, is enough to strongly flatten the BRDF. 


4 Method 


We designed an inversion method aiming at massive data 
analysis. This method consists in two steps : first the genera¬ 
tion of a synthetic database that is representative of the vari¬ 
ability of the model, and then comparison with actual data. 
To generate the synthetic database, we used optical constants 
for water ice at 270 K. The 7K difference between the ac¬ 
tual temperature of the room and the temperature assumed 
for the optical constants has negligible effect. We combined 
the datasets of | Warren and Brandt] ( |2008| ) and |Schmitt et al.| 
( |1998| ) making the junction at 1 /im, the former set for the 
shorter wavelengths, and the latter for the wavelengths larger 
than 1 /im. 

In order to validate the model on the specular reflection 
from the slab, we chose to use the reflectance at 1.5 /im, 
where the ice is very absorbing. Figures [4| and [5] clearly 
demonstrate that there is negligible diffuse contribution in 
geometry outside the specular lobe from the sample with 
12.51mm thick pure slab. Thus, the roughness parameter 6 
is the only one impacting the reflectance in the model. We 
chose to inverse this parameter first, and validate the specu¬ 
lar contribution. 

We will then focus on the validation in the spectral do¬ 
main, for the diffuse contribution. We will use the estimation 
of the roughness parameter 0 obtained earlier and the spec¬ 
tral data in order to estimate the slab thickness and the grain- 
size of the snow substrate. To do this, we assumed that the 
roughness was not changing significantly enough to have a 
notable impact on diffuse reflectance from one sample to an¬ 
other. This assumption is justified by the fact that the differ¬ 
ent columnar ice samples were made the same way, as flat as 
possible and the low value of 0 retrieved as discussed in the 
next section. It is confirmed by the results of section 4.2 that 
suggest a very low roughness, as expected. Such low rough¬ 
ness parameter have negligible influence on the amount of 
energy injected into the surface. 
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4.1 Inversion strategy 

The inversion consists in estimating the model parameter m 
from the model F(m ), that are close to the data d. Taran- 


|tola| ( |1982| ) showed that it can be mathematically solved by 
considering each element as a Probability Density Function 
(PDF). In non-linear direct problem, the solution may not be 
analytically approached. Nevertheless, it is possible to sam¬ 
ple the solution PDF with a Monte Carlo approach as shown 
in fMosegaard and Tarantola (1995), but this solution is very 
time consuming. 

The actual observation is considered as prior information 
on data po{d) in the observation space D. It is assumed to be 
a iV-dimension gaussian PDF Q(d mes ,C), with mean d mes 
and covariance matrix C. The values r % are the observations 
for each element (angular or spectral as described later). The 
covariance matrix C is assumed here to be diagonal since 
measurements at a given geometry/wavelength are indepen¬ 
dent of the other measurements. The diagonal elements Cu 
are o \,... ,cr^, with 07 being the standard deviation. The 
prior information on model parameters pM (m) in the param¬ 
eter space M is independent to the data and corresponds to 
the state of null information if no information is available on 
the parameters. We consider an uniform PDF in their defini¬ 
tion space M. The state of null information represents 

the case when no information is available. It is trivial in our 
case and represent the uniform PDF in the parameters space 
M. The posterior PDF in the model space cjm (m) as defined 
by [Tarantola] ( [ 1982[ ) is : 

°M{m) = kp M (rn)L(m) (1) 

where k is a constant and L(m) is the likelihood function 

, p D {d) 6(d | m) 


L{m) = J 

D 


d d- 


H D (d) 


( 2 ) 


where 0{d\m) is the theoretical relationship of the PDF for 
d given m. We do not consider errors on the model itself, so 
6(d | m) = S(F(m)) also noted for simulated data. So 
the likelihood is simplified into : 


L(m) = G{F(m) - d mes ,C) 


(3) 


and in the case of an uniform prior information pM (to), the 
posterior PDF is: 


(J M (m) = kL(m) (4) 

this expression is explicitly : 

( 1 t =-i 

a M (m) — &.exp ( - - x (F(m) - dmes)C ( F(m ) - dmes) 


The factor k is adjusted to normalize the PDF. The mean 
value of the estimated parameter can be computed by : 


(m) = / 


( 6 ) 


M 


and the standard deviation: 


r (m) 


7 (ra - 

M 


■ m ) 2 .(JM(^)dm 


(7) 


In order to speed up the inversion strategy but keep the ad¬ 
vantage of the Bayesian approach, we choose to sample the 
parameter space M with regular and reasonably fine steps, 
noted i. The likelihood for each element is: 


L(i) — exp ( -- x (d sim(k) dmes) C i^dsimik) dmes) 


( 5 ) 


( 8 ) 


The derivation of posterior PDF with such formalism for 
specular lobe inversion and for spectral inversion are ex¬ 
plained in the next sections. 

4.2 Specular lobe 

To study the specular spot, we have to consider the whole an¬ 
gular sampling of the spot as single data measurement. Simi¬ 
lar to the “pixel” (contraction of picture element ), we choose 
to define the “angel” (contraction of angular element ), as a 
single element in a gridded angular domain. Interestingly, an¬ 
gel also refer to a supernatural being represented in various 
forms of glowing light. A single angel measurement could 
not well constrain the model, even at different wavelengths. 
Instead a full sampling around the specular lobe should be 
enough, even at one single wavelength. We chose a wave¬ 
length where the diffuse contribution was negligible in or¬ 
der to simplify the inversion strategy. We first generated a 
synthetic database (look-up table), using the direct radiative 
transfer model. We simulated spectra in the same geometri¬ 
cal conditions, for a 12.5mm thick ice layer over a granular 
ice substrate constituted of 1000 pm wide grains. These two 
last parameters are not important since the absorption is so 
high in ice, such the main contribution is from the specular 
reflection, and the diffuse contribution is negligible. 

The sampling of the parameters space, that is the lookup 
table, must represent correctly every possible variability. For 
this study, we sampled the roughness parameter from 0 . 1 ° 
to 5° with a constant step dO = 0.01°. We use a likelihood 
function L defined in Eq.[ 8 j where d s i m and d mes are n georn - 
elements vectors, with n georn the number of angel (63 in that 
study). They represent respectively the simulated and mea¬ 
sured reflectance at a given wavelength in every geometry. C 
is a n georn x n georn matrix. It represents the uncertainties on 


°geom 
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the data. In this case, we considered each wavelengths inde¬ 
pendently, thus generating a diagonal matrix, containing the 
level of errors given by the technical data of the instrument 
given by Brissaud et al. (2004). The roughness parameter 0 
returned by the inversion will be described by its normalized 
PDF: 


v{e(i)} 


L (i) d6 


m 

E : LU) 


(9) 


The best match, is the value Q(i) with the highest probability. 
The full PDF can be estimated by its mean : 


<*> 


E Mi) 


and associated standard deviation : 




'E* (*«-<*»■m 

EiL(i) 


( 10 ) 


( 11 ) 


Secondly, we inverted the BRDF as a whole, for each sam¬ 
ple. For this method, and d mes are respectively the 
simulated and measured BRDF, and thus are nb x n geom - 
elements vectors (2379 in that study), where nb is the number 
of bands (61 in that study) and n georn is the number of ge¬ 
ometries (39 in that study) sampled. C is a (n& x n geom ) x 
(rib x Tigeom) diagonal matrix, containing the error on the 
data. We represent the results the same way as previously, but 
there are two parameters to inverse. For the sake of readabil¬ 
ity, we plot the normalized marginal probability density func¬ 
tion for each parameter. We present here the general method 
for the inversion of n p = 2 parameters : the slab thickness 
and the grain-size of the substrate. The PDF for the two pa¬ 
rameters p is described by : 




L(i,j)dp(i,j) 

E i T l j L (hj) d P(h3) 


(13) 


We give error bars on the results that correspond to two stan- For a given parameter p 1 , the marginal PDF of the solution 
dard deviation, and thus a returned value for 0 that is i s ; 


6 r = (0) ±2(J ^ 

4.3 Diffuse spectra 


( 12 ) 


P{pi(i)} 


V (i) dpi (i) 

Ei J2jL(i,j)dp(i,j) 


(14) 


When out of the specular spot, the radiation is controlled by 
the complex transfer through the media (slab ice and bottom 
snow). The experimental samples were made of pure water 
slab ice, without impurity. We generated the look up table 
for every measurement geometry at very high spectral reso¬ 
lution (4. 10“ 2 nm) as required by the variability of the opti¬ 
cal constants of water ice, and then down sampled it at the 
resolution of the instrument (2nm). We sampled the 17,085 
combinations of two parameters for the 39 different geome¬ 
tries : pi the thickness of the slab from 0 mm to 20 mm (noted 
i = [1,201]) every 0.1mm (noted dpi), and p 2 the grain-size 
of the granular substrate from 2 pm to 25 pm every 1 pm and 
from 25 pm to 1500 pm every 25 pm (noted j = [1,85] and 
the corresponding dp 2 (j)). The parameter space is thus ir¬ 
regularly paved with d p(ij) = dpi .dp 2 ( j ). 

For the inversion, we used the same method as previously 
described, with a likelihood function L that is written as 
in equation [8] Two different strategies were adopted. First, 
we inverted each spectra independently. 39 geometries were 
sampled (described in section [T2| ), and thus we conducted 39 
inversions for each sample. This time and d mes are thus 
respectively the simulated and measured spectra. Then 
and d mes are n^-elements vectors, where rib is the number of 
bands (61 in that study) and C is a rib x rib matrix. As pre¬ 
viously (see Section 42), we considered each wavelengths 
independently, thus generating a diagonal matrix, containing 
the level of errors given by the technical data of the instru¬ 
ment given by Brissaud et al. ( |2004| . The error is a percent¬ 


age of the measurement, and thus C will be different for ev¬ 
ery inversion. 


with L' ( i ) = J2j L(i,j)dp 2 (j)- The best match is the value 
Pi(i) with the highest probability. The marginal PDF can be 
described by the mean : 


, v = TnPi (i)L' (i) dpi CO 
W EiE^(U)d P(i,3) 


(15) 


and the associated standard deviation : 


a (pi) ~ 


' E i (Pi (i) ~ (pi)) p_(i) d Pi (i) 

'Ei'Li L (h3)dp(i,j) 


(16) 


As for the roughness parameter, we give error bars on the 
results that correspond to two standard deviation, and thus a 
returned value for p\ that is 


Pir = (pi)±2a {pi) 


(17) 


5 Results 

5.1 Specular lobe 

We performed the inversion taking into account 63 angels 
measurement, but for the sake of readability, figure [4] repre¬ 
sents only the reflectance in the principle plane. The shapes 
and the intensities on Figure [4^ are compatible, but mea¬ 
surement and simulation are not centered at the same point. 
The simulation is centered at the geometrical optics specular 
point (emergence 50° and azimuth 180°), whereas the mea¬ 
surement seems centered around an emergence of 50.5°. This 
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Specular lobe in the principal plane 



emergence angle (°) 


(a) 

Specular lobe in the principal plane 



Figure 4. (a) Measured (black) and simulated (red) reflectance at 
1.5 jum in the principal plan for an incidence angle of 50°. The 
specular lobe measured is not centered at 50°. The sample may be 
slightly misadjusted resulting in a general drift on the observation, 
(b) Measured (black) and simulated (red) reflectance at 1.5 /tm in 
the principal plan for an incidence angle of 50°. We simulated a 
small miss-adjustment of the sample, resulting in a shift of the ob¬ 
servation of 0.5° in emergence and 0.2° in azimuth. With this ad¬ 
justments, the model reproduce well the data. 


could be due to slightly mis-adjustment of the rotation axis 
of the sample in the instrument. This kind of mis-adjustment 
are common, and can easily result in a notable shift up to 1° 
of the measure. We simulated different possible shifts in this 
range, and found a best match represented on Figure [4]} for 
a shift of 0.5° in emergence, as it was suggested by the first 


Measured reflectance at >1=1.5 /jm 



0.0 5.0 10.0 15.0 20.0 25.0 30.0 
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Simulated reflectance at A=1.5 gm 
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Figure 5. Measured and simulated reflectance around the specular 
geometry at 1.5/xm for an incidence angle of 50°. The simulation 
was computed assuming the determined shift of 0.5° in emergence 
and 0.2° in azimuth. The shape as well as the intensity of the spec¬ 
ular lobe are well reproduced. 
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Figure 6. Probability density function a posteriori for the roughness 
parameter 6 , noted V {$}. This function is very sharp and thus the 
parameter 6 is well constrained. The inverted value at 2a is 6 = 
0.424° =b 0.046°. The best match is obtained for 0 — 0.43°. 


plot on figure]?^, and 0.2° in azimuth. The measurements 
and the best match are represented on figure [5] The shape as 
well as the magnitude of the specular lobe are very well re¬ 
produced. Both lobes show a little asymmetry forward. This 
asymmetry is not due to the sampling as it is present also 
when the simulation is not shifted (see the red curve on fig¬ 
ure 0. It is due to an increase of the Fresnel reflexion co¬ 
efficient when the phase angle increase for this range of ge¬ 
ometries. Figure [ 6 ] shows the PDF a posteriori for the pa¬ 
rameter 6. The best match was obtained with 0 = 0.43°. The 
inversion method gives a results close to gaussian shape at 
6 = 0.424° =b 0.046°. Unfortunately, we have no direct mea¬ 
surements of 0. It would require a digital terrain model of 
the sample that is difficult to obtain in icy samples. Still we 
find a low value, that is consistent with the production in lab¬ 
oratory of slabs of columnar ice, that is very flat, but still 
imperfects as described in the dataset. The average slope is 
compatible with a long wavelength slope at the scale of the 
sample, demonstrating that the micro-scale was not impor¬ 
tant in our case. Indeed, for a sample that has a length L , an 
la standard deviation on the thickness Ah can be attributed 
to a general slope i9 = arc tan (^) due to a small error in 
the parallelism of the two surfaces of the slab. In the case of 
sample 3, L = 20 cm and Ah = 1.35mm result in d = 0.39°, 
which is compatible with the roughness given by the inver¬ 
sion. We thus think that what we see is an apparent rough¬ 
ness due to a small general slope on the samples, and that the 
roughness at the surface is much lower than this value. 

Moreover, the value retrieved by the inversion is very well 
constrained as the probability density function is very sharp. 


This means that we have an a posteriori uncertainty on the 
result that is very low. The quality of the reproduction of the 
specular spot by the model suggest that the surface slope de¬ 
scription is a robust description despite its apparent simplic¬ 
ity. Especially, one single slope parameter is enough to de¬ 
scribe this surface. 

5.2 Diffuse 

To reproduce diffuse reflectance we used the results obtained 
with the specular measurement and assumed that the rough¬ 
ness of the samples was not changing much between the ex¬ 
periments. The range of variations in roughness should be 
negligible in the spectral analysis. We simulated slabs over 
snow letting as free parameters the grain-size of the substrate 
and the thickness of the slab. Figure [ 7 ] represent examples of 
the best matches we obtained for the three measured samples 
at various geometries. We also represented the mismatch be¬ 
tween the best fit and the observations. We find an agreement 
between the data and the model that is acceptable. Never¬ 
theless, there seems to be a decrease of quality in the fits as 
the thickness increase. Figure [ 8 ] shows on the same plot an 
example of the marginal PDF for the three samples that are 
associated with the previous fits. Figure [ 8 ] for the thickness 
of the slab (a) and the grain-size of the snow substrate (b). 
The thickness is well constrained as the marginal probability 
density functions a posteriori are relatively sharp and very 
close to gaussian. On the contrary, the grain-size of the sub¬ 
strate seems to have a limited impact on the result since it is 
little constrained. The marginal PDF for the grain-size of the 
substrate are broad, and thus the a posteriori relative uncer¬ 
tainty on the result is very high. Unfortunately, we have no 
reliable measurement of the grain-size of the substrate, as it 
is evolving during the time of the measurements. The general 
trend of decreasing grainsize seems to be in agreement with 
eye’s impression. 

Figure [9] show the measurements and the final result of 
the inversion of the thickness for the three samples, and for 
39 measurement geometries independently. The data and the 
model are compatible. Still, the thickness of the sample 1 is 
slightly overestimated. This may reveal a sensitivity limit of 
the model. The thickness of sample 3 seems underestimated. 
This could be partly due to the duration of the measurement: 
the slab sublimates as the measure is being taken. Moreover, 
the specular measurements were performed on that sample 
increasing even more the duration of the experiment. The 
inversions points on Figure [9] are sorted by increasing inci¬ 
dence, and for each incidence, by increasing azimuth. There 
seems to have an influence of the geometry on the returned 
result : it is particularly clear for sample 2. The estimated 
thickness tend to increase with incidence and decrease with 
azimuth. This effect disappear for large thicknesses (sam¬ 
ple 3). Figure [TO] show the measure and the best match at 
the A = 1.0 pm wavelength, when conducting the inversion 
on the whole BRDF dataset for each sample. The relatively 
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flat behavior of the radiation with the phase angle is reason¬ 
ably well reproduced. The quality of the geometrical match 
increases with the thickness of the sample. This is consis¬ 
tent with the fact that a thicker slab will permit a stronger 
isotropization of the radiation. It is also consistent with the 
disappearance of the geometrical dependence on the estima¬ 
tion for large thicknesses noted on Figure [9] The values of 
thicknesses returned by the inversion are displayed on Fig¬ 
ure EH* : they are also compatible with the data, and the re¬ 
sults are close to the one given by independent inversions on 
each geometry (see Figures[8]and[9]). The grainsizes returned, 
even if compatible with the independent inversion results, are 
at the boundary of the definition range of the parameter, for 
samples 2 and 3. This mean that the model cannot estimate 
this parameter correctly. Indeed, as displayed on Figure |TT] 
the a posteriori marginal probability density functions for the 
samples 2 and 3 are very close to a Dirac at the lower limit of 
the domain. This mean that the model inversion process can¬ 
not fit a value for this parameter inside the definition domain 
that is fully satisfying. This suggest an evolution of the con¬ 
ditions of the experiment between the measure for the sam¬ 
ple 1 and the others. The fact that the returned value is at the 
lower boundary of the grain-size range suggest the the actual 
grainsize of the snow is lower than this value. Unfortunately, 
such grain-size would contradict the fundamental hypothe¬ 
sis of geometrical optics assumed by the model. This results 
thus shall be interpreted as grain-sizes smaller than the limit 
of detection. This kind of very small grain-size could be pro¬ 
duced during the experiments, by a small temperature differ¬ 
ence between the slabs and the natural snow, resulting in the 
condensation of frost at the bottom of the slab layer. 


6 Discussion and conclusion 


The aim of this present work is to validate an approximate 
radiative transfer model developed in Andrie u et al.| ( |2015| ) 
using several assumptions. The most debating one is the 
isotropization of the radiation when it reaches the substrate. 
We first validated qualitatively this assumption with snow 
and ice data. We then quantitatively tested and validated our 
method using a pure slab ice with various thickness and snow 
as a bottom condition. The thicknesses retrieved by the inver¬ 
sion are compatible with the measurements for every geom¬ 
etry, demonstrating the robustness of this method to retrieve 
the slab thickness from spectroscopy only. The result given 
by the inversion of the whole dataset is also compatible with 
the measurements. 

We also validate the angular response of such slab in the 
specular lobe. Unfortunately, it was not possible to measure 
the micro-topography in detail to compare with the retrieved 
data. Nevertheless, we found a very good agreement between 
the simulation and the data. The average slope is compati¬ 
ble with a long wavelength slope at the scale of the sample, 
demonstrating that the micro-scale was not important in our 


case. This is probably due to the sharp slicing method used. 
In future work, an experimental validation of the specular 
lobe and roughness should be addressed. 

The large uncertainties on the grain size inversion demon¬ 
strate that the bottom condition is less important than the slab 
for the radiation field at first order, as expected. Even at a 
thickness of 1.4 mm, since water ice is highly absorbent, the 
bottom layer is difficult to sense. 

On figure [ 7 ] there seems to be a decrease in the quality 
of the fit when the slab thickness increase. We explain it by 
the order in which the experiments were conducted. Indeed, 
the first measurement was of the thinnest slab and the last on 
the thickest. In the meantime, the snow substrate was subli¬ 
mating. The errors in the fits could be due to the increasing 
contamination in the substrate. The natural snow cannot be 
perfectly pure : as it sublimates during the measurements, 
the contribution of the contaminants become stronger and 
stronger. These contaminant are not known, and not taken 
into account in the model. A way to avoid this problem could 
be to set the slab ice on top of a non volatile granular ma¬ 
terial, such as dry mineral sand, which optical constants are 
known or can be determined. However this would not solve 
another problem that is the re-condensation of water into 
frost between the granular substrate and the slab. 

The comparison of the a posteriori uncertainties on the 
thickness of the slab and the grain-size of the snow substrate 
illustrate the fact that those uncertainties depend both on the 
constraint brought by the model itself and the uncertainty 
introduced in the measurement, that only the bayesian ap¬ 
proach can handle. The use of bayesian formalism, is thus 
very powerful in comparison with traditional minimization 
techniques. 

We propose here an fast and innovative inversion method 
aiming at massive inversion, for instance for remote sensing 
spectro-imaging data, that enables to accurately estimate the 
uncertainties on the model’s parameters. As an example, the 
look-up table used for this project were computed in ~ 150 s 
for the roughness study (1763 wavelengths sampled, 30933 
spectra), and ~ 2.5h for the thickness and grainsize study 
(33,186 wavelengths sampled, 666,315 spectra). The inver¬ 
sion itself were performed in less than one tenth of a sec¬ 
ond for specular lobe and independent spectral inversions, 
and 2 s for BRDF-as-a-whole inversions. Every calculation 
was computed on one 4 GB RAM Intel CPU. It has to be 
noted that once the lookup table has been created, an unlim¬ 
ited number of inversion can be conducted. The model is fast 
and the inversion is highly parallel and thus adapted to the 
study of the compact ice covered surfaces of the Solar sys¬ 
tem. For inversion over very large databases, the code has 
been adapted to GPU paralellization. It is also possible to in¬ 
crease de speed of the calculation of the look-up tables by 
doing a multi-CPU computing. 
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Figure 7. Measured and best match of simulated reflectance spectra 
for the geometry of the best match for each sample : at incidence 
40°, emergence 10° and azimuth 140° for sample 1 (a), at incidence 
40°, emergence 20° and azimuth 45° for sample 2 (b), and at inci¬ 
dence 60°, emergence 0° for sample 3 (c). The thicknesses indicated 
were measured before setting the sample in the spectro-goniometer, 
and the errors are given at 2cr. The absolute differences are repre¬ 
sented in blue on each graph. The simulated spectra well reproduce 
the data within the range of a priori uncertainties. For the Sample 3 
(c), the reflectance in the 0.8 pm — 1.0pm range are not very well 
reproduced. The model cannot match the high levels of the mea¬ 
surement. This could be explained by a change in the experimental 
protocol, leading to the condensation of very fine frost at the bottom 
of the slab layer. 
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Figure 8. Marginal probability density functions a posteriori for 
(a) the thickness of the slab V{pi(i)} and (b) the grain-size of the 
snow substrate V {p 2 (j )} for the three samples, and for the geome¬ 
tries described in Figure [ 7 ] The functions are very sharp and very 
close to gaussian for the thickness of the slab (a), but are broad for 
the grain-size of the substrate (b). The thickness is well constrained 
by the inversion, whereas the grain-size of the substrate cannot be 
determined with high precision. 



Figure 9. Results of the inversion and measures with error bars at 
2 a for samples 1, 2 and 3, and for the 39 different geometries of 
measurement. The thicknesses retrieved and measured are compat¬ 
ible. The inversion points (in red) are sorted by incidence (3 val¬ 
ues), and each incidence is then sorted by azimuth (13 values : 1 
for emergence 0° and 6 for each 10° and 20° emergences). The ge¬ 
ometry seems to have an impact on the result for sample 1 and 2. 
The thickness estimated seems to increase with incidence, and de¬ 
crease with azimuth. The geometrical effect seems to disappear for 
big thicknesses. 
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(b) Sample 2. Thickness : 7.45 ± 0.84mm 
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(c) Sample 3. Thickness : 12.51 d= 2.7 mm 

Figure 10. Measured and simulated reflectance factor at A = 1 pm 
for (a) sample 1, (b) sample 2 and (c) sample 3. The simulation 
reproduces, if not perfectly, reasonably the geometrical behavior of 
the surfaces. The quality of the geometrical simulation seems to 
increase with the thickness of the slab. This is consistent with the 
isotropization effect of a slab, that will increase with the thickness. 


■-■-■---r 




L + Measure 
| □ Simulation 

I_._._._I_._._._L 


14 
ft 12 

g 10 


cd 

S 4 

2 
0 

0 5 10 15 20 

Thickness of the slab (mm) 



- Sample 1 : 


Sample 2 : 

1 



- Sample 3 ; 

J- 

— . _ . _ 1 _ . _ 

J _ 1 _ 1 _ 

_ . _ . _ i _ . _ . _ . _ . _ ! 


(a) 


Ph 

Q 

Ph 

13 

a 

' l-H 

W) 

5-i 

cd 


0.04 

c —■—■—■—■—i—■—■— 

- Sample 1 j 



Sample 2 : 

0.03 

F- 

- Sample 3 _ 

0.02 

r b 


0.01 

[ 


0.00 

[ . | 



0 500 1000 1500 

Grainsize of the snow (/xm) 

(b) 

Figure 11. Marginal probability density functions a posteriori for 
(a) the thickness of the slab V {pi (i)} and (b) the grain-size of the 
snow substrate V {p 2 (j)} for the three sample thickness. The func¬ 
tions are very sharp and very close to gaussian for the thickness of 
the slab (a). The a posteriori uncertainties on the results are much 
smaller than the previous ones, because the dataset is bigger and 
thus more constraining. Still, these uncertainties are not fully reli¬ 
able, as the model cannot reproduce perfectly within the a priori 
uncertainties the BRDF (see Figure [Toj. (b) The grainsize can be 
determined or the sample 1, and is consistent with the results on 
inversions of single spectra (see Figure ??). On the contrary, they 
cannot be inverted for sample 2 and 3, as the returned probability 
density function is close to a Dirac at the boundary of the definition 
range. 



















